library(xlsx)
library(plyr)
library(dplyr)
library(tidyr)
library(ARTool)
library(agricolae)
library(ggplot2)
library(ggsignif)
library(ggpubr)
library(gridExtra)

setwd("H://OneDrive//研究//3-小实验TU//5-论文//TU对本外竞争的影响//数据与分析")
                                                                           #设置数据读取与输出文件夹

##################
###  数据输入  ###
##################
rm(list=ls())
grow <- read.xlsx("20220928.xlsx",sheetName = "grow", encoding = "UTF-8")  #生长数据
grow$Treat <- as.factor(grow$Treat)
grow$T <- as.factor(grow$T)
grow$U <- as.factor(grow$U)
grow$C <- as.factor(grow$C)
grow$S <- as.factor(grow$S)

leaf <- read.xlsx("20220928.xlsx",sheetName = "leaf", encoding = "UTF-8")  #叶片数据
leaf$T <- as.factor(leaf$T)
leaf$U <- as.factor(leaf$U)
leaf$C <- as.factor(leaf$C)
leaf$S <- as.factor(leaf$S)

ps <- read.xlsx("20220928.xlsx",sheetName = "ps", encoding = "UTF-8")     #光合数据
ps$T <- as.factor(ps$T)
ps$U <- as.factor(ps$U)
ps$C <- as.factor(ps$C)
ps$S <- as.factor(ps$S)

cnp <- read.xlsx("20220928.xlsx",sheetName = "cnp", encoding = "UTF-8")   #元素数据
cnp$T <- as.factor(cnp$T)
cnp$U <- as.factor(cnp$U)
cnp$C <- as.factor(cnp$C)
cnp$S <- as.factor(cnp$S)

T_treat <- read.xlsx("temperature.xlsx",sheetName = "Sheet1", encoding = "UTF-8")   #元素数据
##################
###  定义函数  ###
##################
interactive.yb <- function(df, Col_DV, Col_IV1, Col_IV2, Col_group) {
  library(plyr)
  
  data_CK <- df[df[,Col_IV1] == "0" & df[,Col_IV2] == "0", c(Col_group, Col_DV)]
  data_1 <- df[df[,Col_IV1] == "1" & df[,Col_IV2] == "0", c(Col_group, Col_DV)]
  data_2 <- df[df[,Col_IV1] == "0" & df[,Col_IV2] == "1", c(Col_group, Col_DV)]
  data_C <- df[df[,Col_IV1] == "1" & df[,Col_IV2] == "1", c(Col_group, Col_DV)]
  data_I <- df[df[,Col_IV1] == "1" & df[,Col_IV2] == "1", c(Col_group, Col_DV)]
  for (i in 1:length(Col_DV)) {
    data_I[,Col_DV[i]] <- (data_1[,Col_DV[i]]+data_2[,Col_DV[i]]-2*data_CK[,Col_DV[i]])/abs(data_CK[,Col_DV[i]])
  }
  out <- data.frame()
  for (i in 1:length(Col_DV)) {
    datac <- ddply(data_I,  c(Col_group), .drop=T,
                   .fun = function(xx, col) {
                     c(N    = length (xx[[col]]),
                       mean = mean   (xx[[col]]),
                       sd   = sd     (xx[[col]])
                     )
                   },
                   Col_DV[i]
    )
    datac$DV   <- Col_DV[i]  
    datac$se   <- datac$sd / sqrt(datac$N) 
    datac$down <- datac$mean - datac$se
    datac$up   <- datac$mean + datac$se
    out <- rbind(out, datac)
  }
  out  <- out[,c("C","S","DV","N","mean","sd","se","down","up")]
  out$label[out$down > 0] <- "+"
  out$label[out$up < 0] <- "-"
  out$label[out$down <= 0 & out$up >= 0] <- "AD"
  return(out)
}

orderPvalue.yb <- function(treatment, means, pvalue, alpha = 0.05) {
  n <- length(means)
  z <- data.frame(treatment, means)
  w <- z[order(z[, 2], decreasing = TRUE), ]
  Order <- as.numeric(rownames(w))
  M <- rep("", n)
  g <-list()
  for (i in 1:n) {
    g[[i]] <- Order
    g1 <- which(pvalue[Order[i], ] < alpha)
    for (j in g1) {
      g[[i]] <- g[[i]][-which(g[[i]] == j)]
    }
  }
  for (i in 1:n) {
    d <- c()
    for (j in 1:length(g[[i]])) {
      for (k in j:length(g[[i]])) {
        S <- pvalue[g[[i]][j],g[[i]][k]] < 0.05
        if(S) {
          d <- unique(c(d, g[[i]][k]))
        }
      }
    }
    for (j in d) {
      g[[i]] <- g[[i]][-which(g[[i]] == j)]
    }
  }
  g <- unique(g)
  for (i in 1:length(g)) {
    for (j in 1:length(g[[i]])) {
      M[g[[i]][j]] <- paste(M[g[[i]][j]], letters[i], sep = "")
    }
  }
  out <- data.frame(treatment, means, groups = M)
  return(out)
}

summary.yb <- function(df, Col_DV, Col_IV, Col_group = NULL, 
                       alpha = 0.05, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
  library(plyr)
  library(dplyr)
  library(ARTool)
  library(agricolae)
  
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = ",", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  result <- list()
  
  for (i in c(1:length(group))) {
    result[[i]] <- list()
    df1 <- df[df$group == group[i],]
    
    for (j in c(1:length(Col_DV))) {
      result[[i]][[j]] <- list()
      model <- paste(Col_DV[j], paste(Col_IV, collapse = "*"), sep = "~")
      length2 <- function (x, na.rm=T) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
      }
      datac <- ddply(df1, Col_IV, .drop=T,
                     .fun = function(xx, col) {
                       c(N    = length2 (xx[[col]], na.rm=T),
                         mean = mean    (xx[[col]], na.rm=T),
                         sd   = sd      (xx[[col]], na.rm=T),
                         min  = min     (xx[[col]], na.rm=T),
                         max  = max     (xx[[col]], na.rm=T)
                       )
                     },
                     Col_DV[j]
      )
      datac$se <- datac$sd / sqrt(datac$N) 
      datac$down <- datac$mean - datac$se
      datac$up <- datac$mean + datac$se
      result[[i]][[j]][[1]] <- datac
      result[[i]][[j]][[1]]$DV <- Col_DV[j]
      result[[i]][[j]][[1]]$group <- group[i]
      
      m <- art(eval(parse(text = paste(Col_DV[j], paste(Col_IV, collapse = "*"), sep = "~"))), data = df1)
      #根据需要进行多因素分析
      result[[i]][[j]][[2]] <- anova(m)                                     #查看交互分析结果
      result[[i]][[j]][[3]] <- art.con(m, eval(parse(text = paste(paste("~ "), paste(Col_IV, collapse = "*")))))%>%
        summary() %>%
        mutate(sig = symnum(p.value, corr = FALSE, na = FALSE,
                            cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                            symbols = c("***", "**", "*", ".", " ")
        )
        )                                                       #查看各单元之间交互作用的具体结果
      
      result[[i]][[j]][[1]] <- unite(result[[i]][[j]][[1]], "IV", all_of(Col_IV), sep = ",", remove = F)
      Q <- matrix(1, ncol = length(result[[i]][[j]][[1]]$IV), nrow = length(result[[i]][[j]][[1]]$IV))
      p <- result[[i]][[j]][[3]]$p.value
      m <- 0
      for (k in 1:(length(result[[i]][[j]][[1]]$IV) - 1)) {
        for (l in (k + 1):length(result[[i]][[j]][[1]]$IV)) {
          m <- m + 1
          Q[k, l] <- p[m]
          Q[l, k] <- p[m]
        }
      }
      groups <- orderPvalue.yb(result[[i]][[j]][[1]]$IV, result[[i]][[j]][[1]]$mean, Q)
      colnames(groups) <- c("IV", "means", "groups")
      result[[i]][[j]][[1]] <- merge(result[[i]][[j]][[1]], groups, by = "IV")
      result[[i]][[j]][[1]] <- result[[i]][[j]][[1]][, c("group", "DV", Col_IV, "N", "mean", "sd", "min", "max", "se", "down", "up", "groups")]
    }
  }
  return(result)
}   #输入: df-数据框，Col_IV-自变量(向量), Col_DV-因变量(可多个, 向量);
    #输出: 列表, 一级-因变量(输入顺序), 二级1-统计描述, 二级2-非参数方差分析F, 二级3-非参数多重比较

transpose.yb <- function(list, S, DV) {
  New <- list()
  for (i in c(1:length(list))) {
    for (j in c(1:length(list[[i]]))) {
      for (k in c(1:length(list[[i]][[j]]))) {
        dataf <- cbind(data.frame(DV = rep(DV[j], length(list[[i]][[j]][[k]][,1]))),
                       data.frame(S  = rep(S[i], length(list[[i]][[j]][[k]][,1]))),
                       list[[i]][[j]][[k]])
        if (i == 1 & j == 1) {
          New[[k]] <- dataf
        } else {
          New[[k]] <- rbind(New[[k]], dataf)
        }
      }
    }
  }
  return(New)
}     #将多重列表转化为一层列表，用于输出excel

fig.yb <- function(list, y, is_C = TRUE) {
  library(ggplot2)
  fig_data <- list()
  for (i in c(1:length(list))) {
    for (j in c(1:length(list[[i]]))) {
      if (i == 1) {
        fig_data[[j]] <- list[[i]][[j]][[1]]
      } else {
        fig_data[[j]] <- rbind(fig_data[[j]], list[[i]][[j]][[1]])
      }
    }
  }
  fig_index <- list()
  if (is_C) {
    for (i in c(1:length(fig_data))) {
      fig_data[[i]]$Treat <- paste(fig_data[[i]]$T, fig_data[[i]]$U, sep = "-")
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "0-0"] <- "CK"
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "1-0"] <- "T"
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "0-1"] <- "U"
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "1-1"] <- "T*U"
      fig_data[[i]]$Treat <- factor(fig_data[[i]]$Treat, levels = c("CK", "T", "U", "T*U"))
      fig_data[[i]]$C <- factor(fig_data[[i]]$C, levels = c("Mix", "Mono"))
      fig_data[[i]]$group <- factor(fig_data[[i]]$group, levels = c("SC", "AA"))
      fig_index[[i]] <- ggplot(fig_data[[i]])+
        geom_pointrange(aes(x = Treat, y = mean, ymin = down, ymax = up, color = C),
                        position = position_dodge(width = -0.5), size = .7)+
        geom_text(aes(x = Treat, y = up, label = groups, group = C),
                  position = position_dodge(width = -0.5), vjust = -.5, size = 5)+
        ylim(min(fig_data[[i]]$down), max(fig_data[[i]]$up)*1.12)+
        facet_grid(. ~ group)+
        theme_bw()+
        labs(x = NULL, y = y[i])+
        theme(legend.position = "none",
              strip.text = element_blank(),
              axis.text.x = element_blank(),
              axis.title.y = element_text(size = 10))
      
    }
  } else{
    for (i in c(1:length(fig_data))) {
      fig_data[[i]]$Treat <- paste(fig_data[[i]]$T, fig_data[[i]]$U, sep = "-")
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "0-0"] <- "CK"
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "1-0"] <- "T"
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "0-1"] <- "U"
      fig_data[[i]]$Treat[fig_data[[i]]$Treat == "1-1"] <- "T*U"
      fig_data[[i]]$Treat <- factor(fig_data[[i]]$Treat, levels = c("CK", "T", "U", "T*U"))
      fig_data[[i]]$S <- factor(fig_data[[i]]$S, levels = c("AA", "SC"))
      fig_index[[i]] <- ggplot(fig_data[[i]])+
        geom_pointrange(aes(x = Treat, y = mean, ymin = down, ymax = up, color = S),
                        position = position_dodge(width = -0.5), size = .7)+
        geom_text(aes(x = Treat, y = up, label = groups, group = S),
                  position = position_dodge(width = -0.5), vjust = -.5, size = 5)+
        ylim(min(fig_data[[i]]$down), max(fig_data[[i]]$up)*1.12)+
        theme_bw()+
        labs(x = NULL, y = y[i])+
        theme(legend.position = "none",
              strip.text = element_blank(),
              axis.text.x = element_blank(),
              axis.title.y = element_text(size = 10))
    }
  }
  return(fig_index)
}

##################
###  数据分析  ###
##################
#二级数据计算
grow$TB <- grow$AB + grow$UB                                         #总生物量
grow$RS <- grow$UB / grow$AB                                         #根冠比
RDI <- merge(ddply(grow[grow$C == "Mix",], c("T", "U", "S", "re"), summarise, B1 = sum(TB)),
             ddply(grow[grow$C == "Mix",], c("T", "U", "re"), summarise, B2 = sum(TB)), by = c("T", "U", "re"))
RDI$RDI <- RDI$B1/RDI$B2
RII <- merge(ddply(grow[grow$C == "Mono",], c("T", "U", "S", "re"), summarise, Bmono = sum(TB)/2),
             ddply(grow[grow$C == "Mix",], c("T", "U", "S", "re"), summarise, Bmix = sum(TB)), 
             by = c("T", "U", "S", "re"))
RII$RII <- (RII$Bmono - RII$Bmix)/(RII$Bmono + RII$Bmix)
compete <- merge(RDI, RII, by = c("T", "U", "S", "re"))[,c("T", "U", "S", "re", "RDI", "RII")]

cnp$CN <- cnp$SC / cnp$SN                                              #碳氮比
cnp$CP <- cnp$SC / cnp$SP                                              #碳磷比
cnp$NP <- cnp$SN / cnp$SP                                              #氮磷比


#统计描述与非参数方差分析
result_grow <- summary.yb(grow, c("TB", "RS"), c("T", "U", "C"), "S")
result_leaf <- summary.yb(leaf, c("leaf_a", "leaf_m"), c("T", "U", "C"), "S")
result_ps <- summary.yb(ps, c("A", "gs", "chl", "Fv.Fm"), c("T", "U", "C"), "S")
result_cnp <- summary.yb(cnp, c("CN", "CP", "NP"), c("T", "U", "C"), "S")

interactive_grow <- interactive.yb(grow, c("TB", "RS"), "T", "U", c("T", "U", "C", "S"))
interactive_leaf <- interactive.yb(leaf, c("leaf_a", "leaf_m"), "T", "U", c("T", "U", "C", "S"))
interactive_ps <- interactive.yb(ps, c("A", "gs", "chl", "Fv.Fm"), "T", "U", c("T", "U", "C", "S"))
interactive_cnp <- interactive.yb(cnp, c("CN", "CP", "NP"), "T", "U", c("T", "U", "C", "S"))
interactive <- rbind(interactive_grow, interactive_leaf, interactive_ps, interactive_cnp)

##################
###  结果输出  ###
##################
#表格输出
Table <- c(transpose.yb(result_grow, c("AA", "SC"), c("TB", "RS")),
           transpose.yb(result_leaf, c("AA", "SC"), c("leaf_a", "leaf_m")),
           transpose.yb(result_ps, c("AA", "SC"), c("A", "gs", "chl", "Fv.Fm")),
           transpose.yb(result_cnp, c("SC", "AA"), c("CN", "CP", "NP")))
for (i in c(1,4,7,10)) {
  Table[[i]]$Treat <- paste(Table[[i]]$T, Table[[i]]$U, sep = "-")
  Table[[i]]$Treat[Table[[i]]$Treat == "0-0"] <- "CK"
  Table[[i]]$Treat[Table[[i]]$Treat == "1-0"] <- "T"
  Table[[i]]$Treat[Table[[i]]$Treat == "0-1"] <- "U"
  Table[[i]]$Treat[Table[[i]]$Treat == "1-1"] <- "T*U"
  Table[[i]]$Treat <- factor(Table[[i]]$Treat, levels = c("CK", "T", "U", "T*U"))
}
sheetname <- c("grow_SD", "grow_ANOVA", "grow_MC",
               "leaf_SD", "leaf_ANOVA", "leaf_MC",
               "ps_SD", "ps_ANOVA", "ps_MC",
               "cnp_SD", "cnp_ANOVA", "cnp_MC")
for (i in 1:12) {
  if (i == 1) {write.xlsx(Table[[i]], "result.xlsx",  sheetName = sheetname[i])}
  else {write.xlsx(Table[[i]], "result.xlsx",  sheetName = sheetname[i], append = TRUE)}
}

#绘图
Fig_plant_data <- data.frame("SC" = "S. canadensis", "AA" = "A. argyi")
row.names(Fig_plant_data) <- ""
Fig_plant <- tableGrob(Fig_plant_data, theme = ttheme_minimal(core=list(fg_params=list(fontface="bold.italic"))), cols = NULL)
Fig_plant$widths <- unit(c(0.1, rep(0.45, 2)), "npc")  #绘制处理标注（以表格形式）

Fig_x_data <- data.frame("CK" = "CK", "T" = "T", "U" = "U", "T × U" = "T × U", "", "CK" = "CK", "T" = "T", "U" = "U", "T × U" = "T × U")
row.names(Fig_x_data) <- ""
Fig_x <- tableGrob(Fig_x_data, theme = ttheme_minimal(core=list(base_size = 6)), cols = NULL)
Fig_x$widths <- unit(c(0.1, rep(0.1056, 4), 0.033, rep(0.1056, 4), 0.022), "npc")
#绘制处理标注（以表格形式）

Fig_x_c_data <- data.frame("CK" = "CK", "T" = "T", "U" = "U", "T × U" = "T × U")
row.names(Fig_x_c_data) <- ""
Fig_x_c <- tableGrob(Fig_x_c_data, theme = ttheme_minimal(core=list(fg_params=list(fontface="bold"))), cols = NULL)
Fig_x_c$widths <- unit(c(0.15, rep(0.20, 4), 0.04), "npc")  #绘制处理标注（以表格形式）

Fig_legend_data <- data.frame(
  C = c("Mix", "Mono"),
  x = c(1, 2),
  y = c(0, 0),
  ymin = c(-1, -1),
  ymax = c(1, 1)
)
Fig_legend <- cowplot::get_legend(
  ggplot(Fig_legend_data)+
    geom_col(aes(x = x, y = y, group = C, fill = C))+
    theme_void()+
    scale_fill_discrete(name = "", 
                        breaks = c("Mono", "Mix"), 
                        labels = c("Monoculture", "Mixed culture"))+
    theme(legend.position="top",
          legend.text = element_text(size = 12, face = "bold"))
)  #绘制图例

Fig_legend_c_data <- data.frame(
  S = c("AA", "SC"),
  x = c(1, 2),
  y = c(0, 0),
  ymin = c(-1, -1),
  ymax = c(1, 1)
)
Fig_legend_c <- cowplot::get_legend(
  ggplot(Fig_legend_c_data)+
    geom_col(aes(x = x, y = y, group = S, fill = S))+
    theme_void()+
    scale_fill_discrete(name = "", 
                        breaks = c("SC", "AA"), 
                        labels = c("S. canadensis", "A. argyi"))+
    theme(legend.position="top",
          legend.text = element_text(size = 12, face = "bold.italic"))
)  #绘制图例

Fig_grow <- fig.yb(result_grow, 
                   c("Total biomass (g)", "Root shoot ratio"))
Fig_leaf <- fig.yb(result_leaf, 
                   c(expression("Leaf area (" * mm ^ 2 * ")"), 
                     "Leaf mass (mg)"))
Fig_ps   <- fig.yb(result_ps, 
                   c(expression("Net photosynthetic (" * μmol ~ m ^ -2 ~ s ^ -1 * ")"),
                     expression("Stomatal conductance (" * μmol ~ m ^ -2 ~ s ^ -1  * ")"), 
                     "Chlorophyll content (SPDA)", "Fv/Fm"))
Fig_cnp  <- fig.yb(result_cnp, c("C:N", "C:P", "N:P"))

Fig_interactive <- ggplot(interactive)+
  geom_pointrange(aes(x = DV, y = mean, ymin = down, ymax = up, color = C),
                  position = position_dodge(width = -0.5), size = .7)+
  geom_text(aes(x = DV, y = up, label = label, group = C),
            position = position_dodge(width = -0.5), vjust = -.5, size = 3)+
  scale_x_discrete(limits = c("TB","RS","leaf_a","leaf_m","A","gs","chl","Fv.Fm","CN","CP","NP"), 
                   label = c("Total biomass", "Root shoot ratio","Leaf area","Leaf mass",
                             "Net photosynthetic","Stomatal conductance", 
                             "Chlorophyll content", "Fv/Fm","C:N", "C:P", "N:P"))+
  scale_color_discrete(name = "", 
                       breaks = c("Mono", "Mix"), 
                       labels = c("Monoculture", "Mixed culture"))+
  facet_grid(S ~ .,labeller=labeller(S = c(SC = "S. canadensis", AA = "A. argyi")))+
  ylim(min(interactive$down), max(interactive$up)*1.10)+
  theme_bw()+
  labs(x = NULL, y = "Relative interaction effect")+
  theme(axis.text.x = element_text(angle = 15,vjust = 0.5, size = 11),
        axis.title.y = element_text(size = 12, face="bold"),
        strip.text.y = element_blank(),
        legend.position="top",
        legend.text = element_text(size = 12, face = "bold"))

interactive_sum <- ddply(interactive, c("S", "C", "label"), summarize, N = length(mean))
interactive_sum$label[interactive_sum$label == "-"] <- "Antagonism (-)"
interactive_sum$label[interactive_sum$label == "+"] <- "Synergies (+)"
interactive_sum$label[interactive_sum$label == "AD"] <- "Additive (AD)"
interactive_sum <- ddply(interactive_sum, c("S", "C"), mutate, N2 = cumsum(N))
interactive_sum$C <- factor(interactive_sum$C, levels = c("Mono", "Mix"))

Fig_interactive_sum <- ggplot(interactive_sum) +
  geom_col(aes(x = 2, y = N, fill = factor(label, levels = c("Additive (AD)", "Synergies (+)", "Antagonism (-)"))), width= 1)+
  geom_text(aes(x = 2, y = N2 - N/2, label = N), size = 5)+
  xlim(0, 2.5)+
  facet_grid(S ~ C,labeller=labeller(S = c(SC = "S. canadensis", AA = "A. argyi"),
                                     C = c(Mono = "Monoculture", Mix = "Mixed culture")))+
  coord_polar(theta="y") +
  theme_void() +
  theme(strip.text.y = element_text(size=12, face="bold.italic", angle = -90),
        strip.text.x = element_text(size=12, face="bold"),
        legend.position ="bottom",
        legend.title=element_blank(),
        legend.text = element_text(size = 10, face = "bold"))


Fig1 <- annotate_figure(
  annotate_figure(
    ggarrange(Fig_grow[[1]], Fig_grow[[2]], Fig_leaf[[1]], Fig_leaf[[2]],
              nrow = 4, ncol = 1, align = "hv"),
    bottom = Fig_x),
  bottom = Fig_legend,
  top = Fig_plant)


Fig2 <- annotate_figure(
  annotate_figure(
    ggarrange(Fig_ps[[1]], Fig_ps[[2]], Fig_ps[[3]], Fig_ps[[4]],
              nrow = 4, ncol = 1, align = "hv"),
    bottom = Fig_x),
  bottom = Fig_legend,
  top = Fig_plant)

Fig3 <- annotate_figure(
  annotate_figure(
    ggarrange(Fig_cnp[[1]], Fig_cnp[[2]], Fig_cnp[[3]],
              nrow = 3, ncol = 1, align = "hv"),
    bottom = Fig_x),
  bottom = Fig_legend,
  top = Fig_plant)

Fig4 <- ggarrange(Fig_interactive, Fig_interactive_sum,
                  nrow = 1, ncol = 2, widths  = c(0.6,0.4))

###存图
ggsave(file = "Fig1.svg", plot = Fig1, width = 7, height = 10,  bg = "white")
ggsave(file = "Fig2.svg", plot = Fig2, width = 7, height = 10,  bg = "white")
ggsave(file = "Fig3.svg", plot = Fig3, width = 7, height = 7.5, bg = "white")
ggsave(file = "Fig4.svg", plot = Fig4, width = 10, height = 5, bg = "white")

ggsave(file = "Fig1.tiff", plot = Fig1, width = 7, height = 10,  bg = "white")
ggsave(file = "Fig2.tiff", plot = Fig2, width = 7, height = 10,  bg = "white")
ggsave(file = "Fig3.tiff", plot = Fig3, width = 7, height = 7.5, bg = "white")
ggsave(file = "Fig4.tiff", plot = Fig4, width = 8, height = 5, bg = "white")

ggsave(file = "Fig1.png", plot = Fig1, width = 7, height = 10,  bg = "white")
ggsave(file = "Fig2.png", plot = Fig2, width = 7, height = 10,  bg = "white")
ggsave(file = "Fig3.png", plot = Fig3, width = 7, height = 7.5, bg = "white")
ggsave(file = "Fig4.png", plot = Fig4, width = 10, height = 5, bg = "white")

